############################################################################################# Preamble ----
# File:         nomikos_jop_2021_replication.R
# Description:  Replication script for "Peacekeeping and the Enforcement of Intergroup 
#               Cooperation: Evidence from Mali," Journal of Politics, 2021.
# Author:       William G. Nomikos
# Last Update:  1/29/21
# Instructions: Download the "Replication" folder from the JoP dataverse and run this 
#               .R script to produce all the figures and the tables, which will be 
#               saved in the same folder. 

# **Presets----
rm(list=ls())
setwd("/Users/williamnomikos/Dropbox/Work/Projects/Intergroup Coop/JoP FINAL/Replication") #Replace this with your working directory
getwd()

# **Libraries ----
# If you don't have any of these libraries, make sure to run "install.packages'' first
library(sandwich)
library(arm)
library(readxl)
library(haven)
library(tidyverse)
library(xtable)
library(knitr)
library(kableExtra)
library(multiwayvcov)
library(sandwich)
library(lmtest)
library(stargazer)
library(konfound)
library(plotrix)
library(sciplot)

# **Data Management---- 

nomikos_jop_data<-read_csv("nomikos_jop_data.csv")

# Categorizing uncontact
table(nomikos_jop_data$un_contact)
lab_high_uncontact<-subset(nomikos_jop_data,)
lab_med_uncontact<-subset(nomikos_jop_data,un_see>=1)
lab_low_uncontact<-subset(nomikos_jop_data,un_talk==0)


#By Self-Reported Gender
lab_male<-subset(nomikos_jop_data,female==0)
lab_female<-subset(nomikos_jop_data,female==1)


############################################################################################## Article ----

# **Figure 1: Main Results ----

# ** Full Sample 
# **** Control group 
model<-lm(donation~fra+un+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg<-coeftest(model,vcov)
reg


new.data<-data.frame(group=c("C","UN","Fr"))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# **** UN group 
model<-lm(donation~untreat+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg<-coeftest(model,vcov)
reg

df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)


# **** France group 
model<-lm(donation~fratreat+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg<-coeftest(model,vcov)
reg

df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[3]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(reg[1,1]+est,0)

ggplot(aes(x=group, y=fit,ymin=ymin95, ymax=ymax95),data=new.data)+
  geom_linerange(size=1)+
  theme_bw()+
  ylab("Amount Sent to Tuareg Partner") + xlab("Treatment") +
  stat_summary(fun.data = mean_cl_normal, fun.args = list(conf.int = 0.95), geom = "errorbar",
               ymin=new.data$ymin95,ymax=new.data$ymax95, position = position_dodge(width = 0.90), width = 0.1) +
  geom_point(size=18,shape=21, fill=c("white", "lightblue", "grey"))+
  geom_text(x=c(1,2,3), y=new.data$fit, label=round(new.data$fit,2), size=5)+
  theme(text = element_text(size=14),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
        axis.ticks.x=element_blank(),
        axis.text.x = element_text(size=12))



dev.copy(pdf,"fig_lab_main.pdf", height=6, width=8)
dev.off()


# **Figure 2: Heterogenous Treatment Effects ----


# TuaTrust 1 

# Regressions
model<-lm(donation~fra*tua_trust_1+un*tua_trust_1+adults_mix+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_control<-coeftest(model,vcov)

model<-lm(donation~untreat*tua_trust_1+adults_mix+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_un<-coeftest(model,vcov)

model<-lm(donation~fratreat*tua_trust_1+adults_mix+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_fra<-coeftest(model,vcov)


# Tuareg Trust 


new.data<-data.frame(subgroup=rep("Tuareg Do Not Support Separatists",3))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))


# **** Control group 
reg<-reg_control
df<-model$df.residual
se<-reg[1,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
new.data$ymin95<-control_est+tuatrust_est-qt(c(.975), df=df)*se
new.data$ymax95<-control_est+tuatrust_est+qt(c(.975), df=df)*se
new.data$fit<-round(control_est+tuatrust_est,0)

# **** UN group 
reg<-reg_un
df<-model$df.residual
se<-reg[2,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
un_est<-reg[2,1]
tuatrust_un_est<-reg[19,1]
new.data$ymin95[2]<-control_est+tuatrust_est+un_est+tuatrust_un_est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-control_est+tuatrust_est+un_est+tuatrust_un_est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(control_est+tuatrust_est+un_est+tuatrust_un_est,0)

# **** France group 
reg<-reg_fra
df<-model$df.residual
se<-reg[2,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
fra_est<-reg[2,1]
tuatrust_fra_est<-reg[19,1]
new.data$ymin95[3]<-control_est+tuatrust_est+fra_est+tuatrust_fra_est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-control_est+tuatrust_est+fra_est+tuatrust_fra_est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(control_est+tuatrust_est+fra_est+tuatrust_fra_est,0)


# Tuareg No Trust

old.data<-new.data

new.data<-data.frame(subgroup=rep("Tuareg Support Separatists",3))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))

# **** Control group 

reg<-reg_control
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# **** UN group 

reg<-reg_un
df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)


# **** France group 
reg<-reg_fra
df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[3]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(reg[1,1]+est,0)


# Figure

new.data<-rbind(new.data,old.data)
tua_trust_df<-new.data
# ** TuaTrust 2


# Regressions
model<-lm(donation~fra*tua_trust_2+un*tua_trust_2+adults_mix+children+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_control<-coeftest(model,vcov)

model<-lm(donation~untreat*tua_trust_2+adults_mix+children+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_un<-coeftest(model,vcov)

model<-lm(donation~fratreat*tua_trust_2+adults_mix+children+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_fra<-coeftest(model,vcov)


# Tuareg Trust 


new.data<-data.frame(subgroup=rep("Tuareg are not Biased",3))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))


# **** Control group 
reg<-reg_control
df<-model$df.residual
se<-reg[1,2]
control_est<-reg[1,1]

tuatrust_est<-reg[3,1]
new.data$ymin95<-control_est+tuatrust_est-qt(c(.975), df=df)*se
new.data$ymax95<-control_est+tuatrust_est+qt(c(.975), df=df)*se
new.data$fit<-round(control_est+tuatrust_est,0)

# **** UN group 
reg<-reg_un
df<-model$df.residual
se<-reg[2,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
un_est<-reg[2,1]
tuatrust_un_est<-reg[20,1]
new.data$ymin95[2]<-control_est+tuatrust_est+un_est+tuatrust_un_est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-control_est+tuatrust_est+un_est+tuatrust_un_est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(control_est+tuatrust_est+un_est+tuatrust_un_est,0)

# **** France group 
reg<-reg_fra
df<-model$df.residual
se<-reg[2,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
fra_est<-reg[2,1]
tuatrust_fra_est<-reg[20,1]
new.data$ymin95[3]<-control_est+tuatrust_est+fra_est+tuatrust_fra_est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-control_est+tuatrust_est+fra_est+tuatrust_fra_est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(control_est+tuatrust_est+fra_est+tuatrust_fra_est,0)


# Tuareg No Trust

old.data<-new.data

new.data<-data.frame(subgroup=rep("Tuareg are Biased",3))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))

# **** Control group 

reg<-reg_control
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# **** UN group 

reg<-reg_un
df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)


# **** France group 
reg<-reg_fra
df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[3]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(reg[1,1]+est,0)


# Figure

new.data<-rbind(new.data,old.data)
tua_trust2_df<-new.data





# ** TuaFriend


# Regressions
model<-lm(donation~fra*tua_friend+un*tua_friend+adults_mix+educ+female+neighborhood_makeup+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_control<-coeftest(model,vcov)

model<-lm(donation~untreat*tua_friend+adults_mix+educ+female+neighborhood_makeup+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_un<-coeftest(model,vcov)

model<-lm(donation~fratreat*tua_friend+adults_mix+educ+female+neighborhood_makeup+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_fra<-coeftest(model,vcov)


# Tuareg Trust 


new.data<-data.frame(subgroup=rep("Has Tuareg Friend",3))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))


# **** Control group 
reg<-reg_control
df<-model$df.residual
se<-reg[1,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
new.data$ymin95<-control_est+tuatrust_est-qt(c(.975), df=df)*se
new.data$ymax95<-control_est+tuatrust_est+qt(c(.975), df=df)*se
new.data$fit<-round(control_est+tuatrust_est,0)

# **** UN group 
reg<-reg_un
df<-model$df.residual
se<-reg[2,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
un_est<-reg[2,1]
tuatrust_un_est<-reg[22,1]
new.data$ymin95[2]<-control_est+tuatrust_est+un_est+tuatrust_un_est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-control_est+tuatrust_est+un_est+tuatrust_un_est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(control_est+tuatrust_est+un_est+tuatrust_un_est,0)

# **** France group 
reg<-reg_fra
df<-model$df.residual
se<-reg[2,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
fra_est<-reg[2,1]
tuatrust_fra_est<-reg[22,1]
new.data$ymin95[3]<-control_est+tuatrust_est+fra_est+tuatrust_fra_est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-control_est+tuatrust_est+fra_est+tuatrust_fra_est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(control_est+tuatrust_est+fra_est+tuatrust_fra_est,0)


# Tuareg No Trust

old.data<-new.data

new.data<-data.frame(subgroup=rep("No Tuareg Friend",3))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))

# **** Control group 

reg<-reg_control
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# **** UN group 

reg<-reg_un
df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)


# **** France group 
reg<-reg_fra
df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[3]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(reg[1,1]+est,0)


# Figure

new.data<-rbind(new.data,old.data)
tua_friend_df<-new.data


# ** GovTrust


# Regressions
model<-lm(donation~fra*gov_trust+un*gov_trust+bamako+untrust+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_control<-coeftest(model,vcov)

model<-lm(donation~untreat*gov_trust+bamako+untrust+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_un<-coeftest(model,vcov)

model<-lm(donation~fratreat*gov_trust+untrust+bamako+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg_fra<-coeftest(model,vcov)


# Tuareg Trust 


new.data<-data.frame(subgroup=rep("Trust Government",3))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))


# **** Control group 
reg<-reg_control
df<-model$df.residual
se<-reg[1,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
new.data$ymin95<-control_est+tuatrust_est-qt(c(.975), df=df)*se
new.data$ymax95<-control_est+tuatrust_est+qt(c(.975), df=df)*se
new.data$fit<-round(control_est+tuatrust_est,0)

# **** UN group 
reg<-reg_un
df<-model$df.residual
se<-reg[2,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
un_est<-reg[2,1]
tuatrust_un_est<-reg[19,1]
new.data$ymin95[2]<-control_est+tuatrust_est+un_est+tuatrust_un_est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-control_est+tuatrust_est+un_est+tuatrust_un_est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(control_est+tuatrust_est+un_est+tuatrust_un_est,0)

# **** France group 
reg<-reg_fra
df<-model$df.residual
se<-reg[2,2]
control_est<-reg[1,1]
tuatrust_est<-reg[3,1]
fra_est<-reg[2,1]
tuatrust_fra_est<-reg[19,1]
new.data$ymin95[3]<-control_est+tuatrust_est+fra_est+tuatrust_fra_est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-control_est+tuatrust_est+fra_est+tuatrust_fra_est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(control_est+tuatrust_est+fra_est+tuatrust_fra_est,0)


# Tuareg No Trust

old.data<-new.data

new.data<-data.frame(subgroup=rep("Do Not Trust Government",3))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))

# **** Control group 

reg<-reg_control
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# **** UN group 

reg<-reg_un
df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)


# **** France group 
reg<-reg_fra
df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[3]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(reg[1,1]+est,0)



new.data<-rbind(new.data,old.data)
gov_trust_df<-new.data

# ** Figure

subgroup_treats_df<-rbind(tua_trust_df,tua_trust2_df,tua_friend_df,gov_trust_df)


ggplot(aes(x=group, y=fit,ymin=ymin95, ymax=ymax95),data=subgroup_treats_df,fill=subgroup)+
  scale_color_manual(values=c("#FF3300", "#3399FF", "#33CC00"))+
  geom_bar(stat = "summary", fun.y = "mean", fill = rep(c('grey', 'lightblue', 'darkred'),8)) +
  theme_bw() +
  facet_wrap(~subgroup, ncol = 2)+
  stat_summary(fun.data = mean_cl_normal, fun.args = list(conf.int = 0.95), geom = "errorbar",
               ymin=subgroup_treats_df$ymin95,ymax=subgroup_treats_df$ymax95, 
               position = position_dodge(width = 0.90), width = 0.1) +
  #geom_point(size=12,shape=21, fill="white")+
  geom_text(x=rep(c(1,2,3),8),y=rep(rep(200,3),8),label=rep(c("Control","UN","France"),8))+
  theme(text = element_text(size=14),
        axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0)),
        axis.title.x = element_blank(),
        axis.ticks.x=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank(),
        axis.text.x = element_blank())+
  #geom_point(size=8,shape=22, fill="white")+
  geom_label(label=subgroup_treats_df$fit,size=4,fill="white")+
  ylab("Amount Sent to Tuareg Partner") 



dev.copy(pdf,"fig_lab_het_effects.pdf", height=8, width=6)
dev.off()

# **Figure 3: UN Contact ----

# ** UN see low 

# **** Control group 
model<-lm(donation~fra+un+female+neighborhood_makeup+rebelviol_any+factor(enumid)+factor(day), data=lab_low_uncontact)
summary(model)
vcov <- cluster.vcov(model, lab_low_uncontact$cluster)
reg<-coeftest(model,vcov)
reg

new.data<-data.frame(subgroup=rep("Infrequent Contact",2))
new.data$group<-factor(c("Control", "UN"),levels=c("Control", "UN"))
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# **** UN group 
model<-lm(donation~untreat+female+neighborhood_makeup+rebelviol_any+factor(enumid)+factor(day), data=lab_low_uncontact)
summary(model)
vcov <- cluster.vcov(model, lab_low_uncontact$cluster)
reg<-coeftest(model,vcov)
reg

df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)

new.data1<-new.data

# ** UN see medium
model<-lm(donation~fra+un+female+neighborhood_makeup+rebelviol_any+factor(enumid)+factor(day), data=lab_med_uncontact)
summary(model)
vcov <- cluster.vcov(model, lab_med_uncontact$cluster)
reg<-coeftest(model,vcov)
reg



new.data<-data.frame(subgroup=rep("Some Contact",2))
new.data$group<-factor(c("Control", "UN"),levels=c("Control", "UN"))
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# **** UN group 
model<-lm(donation~untreat+female+neighborhood_makeup+rebelviol_any+factor(enumid)+factor(day), data=lab_med_uncontact)
summary(model)
vcov <- cluster.vcov(model, lab_med_uncontact$cluster)
reg<-coeftest(model,vcov)
reg

df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)

new.data2<-new.data

# ** UN see high

# **** Control group 
model<-lm(donation~fra+un+female+neighborhood_makeup+rebelviol_any+factor(enumid)+factor(day), data=lab_high_uncontact)
summary(model)
vcov <- cluster.vcov(model, lab_high_uncontact$cluster)
reg<-coeftest(model,vcov)
reg

new.data<-data.frame(subgroup=rep("Frequent Contact",2))
new.data$group<-factor(c("Control", "UN"),levels=c("Control", "UN"))
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# **** UN group 
model<-lm(donation~untreat+female+neighborhood_makeup+rebelviol_any+factor(enumid)+factor(day), data=lab_high_uncontact)
summary(model)
vcov <- cluster.vcov(model, lab_high_uncontact$cluster)
reg<-coeftest(model,vcov)
reg

df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)

new.data3<-new.data

#** Draw figure

new.data<-rbind(new.data1,new.data2,new.data3)
new.data$subgroup

effect<-list("30% increase","", "31% increase","","68% increase","")


ggplot(aes(x=group, y=fit,ymin=ymin95, ymax=ymax95),data=new.data)+
  scale_color_manual(values=c("#FF3300", "#3399FF", "#33CC00"))+
  geom_bar(stat = "summary", fun.y = "mean", fill = rep(c('grey', 'lightblue'),3)) +
  theme_minimal() +
  theme_bw()+
  ylab("Amount Sent to Tuareg Partner") + xlab("Treatment") +
  stat_summary(fun.data = mean_cl_normal, fun.args = list(conf.int = 0.95), geom = "errorbar",
               ymin=new.data$ymin95,ymax=new.data$ymax95, 
               position = position_dodge(width = 0.90), width = 0.1) +
  facet_wrap(~subgroup, ncol = 3)+
  geom_point(size=4,shape=21, fill="black")+
  theme(text = element_text(size=14))+
  geom_text(x=rep(1.5,6), y=rep(950,6), label=effect, size=4.5, color="black")+
  ylim(0,1000)


dev.copy(pdf,"fig_lab_un_contact.pdf", height=6, width=10)
dev.off()




############################################################################################## Appendix ----


# **Table B1: Summary Statistics ----
# Create function for easy summary statistics
summ.stats.bal<-function(x,row){
  mean_con <- mean(x[nomikos_jop_data$treatment=="control"], na.rm=T)
  mean_fr <- mean(x[nomikos_jop_data$treatment=="france"], na.rm=T)
  mean_un <- mean(x[nomikos_jop_data$treatment=="un"], na.rm=T)
  
  diff_fr_con <-mean_fr-mean_con
  diff_un_con <-mean_un-mean_con
  diff_fr_un <-mean_fr-mean_un
  
  pval_fr_con<-as.numeric(t.test(x[nomikos_jop_data$treatment=="france"], x[nomikos_jop_data$treatment=="control"], alternative = "two.sided", var.equal = FALSE)[3])
  pval_un_con<-as.numeric(t.test(x[nomikos_jop_data$treatment=="un"], x[nomikos_jop_data$treatment=="control"], alternative = "two.sided", var.equal = FALSE)[3])
  pval_fr_un<-as.numeric(t.test(x[nomikos_jop_data$treatment=="france"], x[nomikos_jop_data$treatment=="un"], alternative = "two.sided", var.equal = FALSE)[3])
  
  t <- round(cbind(mean_con,mean_fr,mean_un,diff_fr_con,diff_un_con,diff_fr_un,pval_fr_con,pval_un_con,pval_fr_un),digits=2)
  tab_balance[row,]<-t
  tab_balance<<-tab_balance
}


tab_balance<-matrix(NA,10,9)
summ.stats.bal(nomikos_jop_data$age,1)
summ.stats.bal(nomikos_jop_data$female,2)
summ.stats.bal(nomikos_jop_data$bambara,3)
summ.stats.bal(nomikos_jop_data$bamako,4)
summ.stats.bal(nomikos_jop_data$educ,5)
summ.stats.bal(nomikos_jop_data$children,6)
summ.stats.bal(nomikos_jop_data$untrust,7)
summ.stats.bal(nomikos_jop_data$adults_mix,8)
summ.stats.bal(nomikos_jop_data$neighborhood_makeup,9)
summ.stats.bal(nomikos_jop_data$rebelviol_any,10)


rownames(tab_balance)<-c("Age","Female","Bambara","Bamako Native","Education (0-8)","Children","UN Trust (0-2)","Adults Mix (0-2)",
                         "Neighborhood Makeup","Exposure to Violence")
colnames(tab_balance) <- c("Control (C)","France (Fr)","UN", 
                           "Fr-C","UN-C","Fr-UN",
                           "Fr-C", "UN-C","Fr-UN")
tab_balance

tab_balance%>%
  kable("latex",booktabs=T,linesep = "", # latex settings
        caption = "Summary Statistics and Balance on Demographic Covariates between Treatments" # caption
  ) %>%
  kable_styling(latex_options = "striped", font_size = 10)%>%
  add_header_above(c(" " = 1, "Mean" = 3, "Difference" = 3, "p-value" = 3)) %>%
  cat(file="tab_summ_stats.tex")






# ** Table D1: Summ Stats Trust (Separatists) ----
balance.viol.opinion<-function(x,row){
  mean_trust <- mean(x[nomikos_jop_data$viol_opinion<=0], na.rm=T)
  mean_notrust <- mean(x[nomikos_jop_data$viol_opinion>=1], na.rm=T)
  
  diff <-mean_trust-mean_notrust
  
  pval<-as.numeric(t.test(x[nomikos_jop_data$viol_opinion<=0], x[nomikos_jop_data$viol_opinion>=1], alternative = "two.sided", var.equal = FALSE)[3])
  
  t <- round(cbind(mean_trust,mean_notrust,diff,pval),digits=2)
  tab_balance[row,]<-t
  tab_balance<<-tab_balance
}


tab_balance<-matrix(NA,10,4)
balance.viol.opinion(nomikos_jop_data$age,1)
balance.viol.opinion(nomikos_jop_data$female,2)
balance.viol.opinion(nomikos_jop_data$bambara,3)
balance.viol.opinion(nomikos_jop_data$bamako,4)
balance.viol.opinion(nomikos_jop_data$educ,5)
balance.viol.opinion(nomikos_jop_data$children,6)
balance.viol.opinion(nomikos_jop_data$untrust,7)
balance.viol.opinion(nomikos_jop_data$adults_mix,8)
balance.viol.opinion(nomikos_jop_data$neighborhood_makeup,9)
balance.viol.opinion(nomikos_jop_data$rebelviol_any,10)

rownames(tab_balance)<-c("Age","Female","Bambara","Bamako Native","Education (0-8)","Children","UN Trust (0-2)","Adults Mix (0-2)","Neighborhood Makeup","Rebelviolence")
colnames(tab_balance) <- c("High Trust","Low Trust", "Difference","p-value")
tab_balance



tab_balance%>%
  kable("latex",booktabs=T,linesep = "", # latex settings
        caption = "Summary Statistics and Balance on Demographic Covariates for Trust (Tuareg Support Separatists)" # caption
  ) %>%
  kable_styling(latex_options = c("striped", "hold_position"), font_size = 10)%>%
  add_header_above(c(" " = 1, "Mean" = 2, "Difference" = 2)) %>%
  cat(file="tab_trust_bal1.tex")



# ** Table D2: Summ Stats Trust (Tuareg friend) ----


balance.tua.friend<-function(x,row){
  mean_trust <- mean(x[nomikos_jop_data$tua_friend==1], na.rm=T)
  mean_notrust <- mean(x[nomikos_jop_data$tua_friend==0], na.rm=T)
  
  diff <-mean_trust-mean_notrust
  
  pval<-as.numeric(t.test(x[nomikos_jop_data$tua_friend==1], x[nomikos_jop_data$tua_friend==0], alternative = "two.sided", var.equal = FALSE)[3])
  
  
  
  t <- round(cbind(mean_trust,mean_notrust,diff,pval),digits=2)
  tab_balance[row,]<-t
  tab_balance<<-tab_balance
}



tab_balance<-matrix(NA,10,4)
balance.tua.friend(nomikos_jop_data$age,1)
balance.tua.friend(nomikos_jop_data$female,2)
balance.tua.friend(nomikos_jop_data$bambara,3)
balance.tua.friend(nomikos_jop_data$bamako,4)
balance.tua.friend(nomikos_jop_data$educ,5)
balance.tua.friend(nomikos_jop_data$children,6)
balance.tua.friend(nomikos_jop_data$untrust,7)
balance.tua.friend(nomikos_jop_data$adults_mix,8)
balance.tua.friend(nomikos_jop_data$neighborhood_makeup,9)
balance.tua.friend(nomikos_jop_data$rebelviol_any,10)


rownames(tab_balance)<-c("Age","Female","Bambara","Bamako Native","Education (0-8)","Children","UN Trust (0-2)","Adults Mix (0-2)","Neighborhood Makeup","Rebelviolence")
colnames(tab_balance) <- c("High Trust","Low Trust", "Difference","p-value")

tab_balance


tab_balance%>%
  kable("latex",booktabs=T,linesep = "", # latex settings
        caption = "Summary Statistics and Balance on Demographic Covariates for Trust (Has Close Tuareg Friend)" # caption
  ) %>%
  kable_styling(latex_options = c("striped", "hold_position"), font_size = 10)%>%
  add_header_above(c(" " = 1, "Mean" = 2, "Difference" = 2)) %>%
  cat(file="tab_trust_bal2.tex")






# ** Table D3: Summ Stats Trust (Tuareg biased) ----


balance.tua.bias.opinion<-function(x,row){
  mean_trust <- mean(x[nomikos_jop_data$bias_opinion==0], na.rm=T)
  mean_notrust <- mean(x[nomikos_jop_data$bias_opinion>=1], na.rm=T)
  
  diff <-mean_trust-mean_notrust
  
  pval<-as.numeric(t.test(x[nomikos_jop_data$bias_opinion==0], x[nomikos_jop_data$bias_opinion>=1], alternative = "two.sided", var.equal = FALSE)[3])
  
  t <- round(cbind(mean_trust,mean_notrust,diff,pval),digits=2)
  tab_balance[row,]<-t
  tab_balance<<-tab_balance
}







tab_balance<-matrix(NA,10,4)
balance.tua.bias.opinion(nomikos_jop_data$age,1)
balance.tua.bias.opinion(nomikos_jop_data$female,2)
balance.tua.bias.opinion(nomikos_jop_data$bambara,3)
balance.tua.bias.opinion(nomikos_jop_data$bamako,4)
balance.tua.bias.opinion(nomikos_jop_data$educ,5)
balance.tua.bias.opinion(nomikos_jop_data$children,6)
balance.tua.bias.opinion(nomikos_jop_data$untrust,7)
balance.tua.bias.opinion(nomikos_jop_data$adults_mix,8)
balance.tua.bias.opinion(nomikos_jop_data$neighborhood_makeup,9)
balance.tua.bias.opinion(nomikos_jop_data$rebelviol_any,10)


rownames(tab_balance)<-c("Age","Female","Bambara","Bamako Native","Education (0-8)","Children","UN Trust (0-2)","Adults Mix (0-2)","Neighborhood Makeup","Rebelviolence")
colnames(tab_balance) <- c("High Trust","Low Trust", "Difference","p-value")
tab_balance



tab_balance%>%
  kable("latex",booktabs=T,linesep = "", # latex settings
        caption = "Summary Statistics and Balance on Demographic Covariates for Trust (Believes Tuareg Are Biased)" # caption
  ) %>%
  kable_styling(latex_options = c("striped", "hold_position"), font_size = 10)%>%
  add_header_above(c(" " = 1, "Mean" = 2, "Difference" = 2)) %>%
  cat(file="tab_trust_bal3.tex")

# ** Table D4: Summ Stats Trust (gov trust) ----

balance.gov.trust<-function(x,row){
  mean_trust <- mean(x[nomikos_jop_data$locgovtrust>=1|nomikos_jop_data$natgovtrust>=1], na.rm=T)
  mean_notrust <- mean(x[nomikos_jop_data$locgovtrust==0|nomikos_jop_data$natgovtrust==0], na.rm=T)
  
  diff <-mean_trust-mean_notrust
  
  pval<-as.numeric(t.test(x[nomikos_jop_data$gov_trust==0], x[nomikos_jop_data$gov_trust>=1], alternative = "two.sided", var.equal = FALSE)[3])
  
  t <- round(cbind(mean_trust,mean_notrust,diff,pval),digits=2)
  tab_balance[row,]<-t
  tab_balance<<-tab_balance
}



tab_balance<-matrix(NA,10,4)
balance.gov.trust(nomikos_jop_data$age,1)
balance.gov.trust(nomikos_jop_data$female,2)
balance.gov.trust(nomikos_jop_data$bambara,3)
balance.gov.trust(nomikos_jop_data$bamako,4)
balance.gov.trust(nomikos_jop_data$educ,5)
balance.gov.trust(nomikos_jop_data$children,6)
balance.gov.trust(nomikos_jop_data$untrust,7)
balance.gov.trust(nomikos_jop_data$adults_mix,8)
balance.gov.trust(nomikos_jop_data$neighborhood_makeup,9)
balance.gov.trust(nomikos_jop_data$rebelviol_any,10)

names(lab_data)

rownames(tab_balance)<-c("Age","Female","Bambara","Bamako Native","Education (0-8)","Children","UN Trust (0-2)","Adults Mix (0-2)","Neighborhood Mix (0-2)","Victimization")
colnames(tab_balance) <- c("High Trust","Low Trust", "Difference","p-value")
tab_balance


tab_balance%>%
  kable("latex",booktabs=T,linesep = "", # latex settings
        caption = "Summary Statistics and Balance on Demographic Covariates for Trust (Trust in Government)" # caption
  ) %>%
  kable_styling(latex_options = c("striped", "hold_position"), font_size = 10)%>%
  add_header_above(c(" " = 1, "Mean" = 2, "Difference" = 2)) %>%
  cat(file="tab_trust_bal4.tex")

# ** Table D4: Summ Stats UN Contact ----


un.contact.bal<-function(x,row){
  mean_low <- mean(x[lab_small$un_contact==0], na.rm=T)
  mean_med <- mean(x[lab_small$un_contact==1], na.rm=T)
  mean_hi <- mean(x[lab_small$un_contact==2], na.rm=T)
  
  diff_med_low <-mean_med-mean_low
  diff_hi_low <-mean_hi-mean_low
  diff_med_hi <-mean_med-mean_hi
  
  pval_med_low<-as.numeric(t.test(x[lab_small$un_contact==1], x[lab_small$un_contact==0], alternative = "two.sided", var.equal = FALSE)[3])
  pval_hi_low<-as.numeric(t.test(x[lab_small$un_contact==2], x[lab_small$un_contact==0], alternative = "two.sided", var.equal = FALSE)[3])
  pval_med_hi<-as.numeric(t.test(x[lab_small$un_contact==1], x[lab_small$un_contact==2], alternative = "two.sided", var.equal = FALSE)[3])
  
  t <- round(cbind(mean_low,mean_med,mean_hi,diff_med_low,diff_hi_low,diff_med_hi,pval_med_low,pval_hi_low,pval_med_hi),digits=2)
  tab_balance[row,]<-t
  tab_balance<<-tab_balance
}




tab_balance<-matrix(NA,10,9)
un.contact.bal(lab_small$age,1)
un.contact.bal(lab_small$female,2)
un.contact.bal(lab_small$bambara,3)
un.contact.bal(lab_small$bamako,4)
un.contact.bal(lab_small$educ,5)
un.contact.bal(lab_small$children,6)
un.contact.bal(lab_small$untrust,7)
un.contact.bal(lab_small$adults_mix,8)
un.contact.bal(lab_small$neighborhood_makeup,9)
un.contact.bal(lab_small$rebelviol_any,10)


rownames(tab_balance)<-c("Age","Female","Bambara","Bamako Native","Education (0-8)","Children","UN Trust (0-2)","Adults Mix (0-2)","Neighborhood Makeup","Rebelviolence")
colnames(tab_balance) <- c("Low","Medium","High", 
                           "Med-Low","High-Low","Medium-High",
                           "Med-Low","High-Low","Medium-High")
tab_balance

tab_balance%>%
  kable("latex",booktabs=T,linesep = "", # latex settings
        caption = "Summary Statistics and Balance on Demographic Covariates for UN Contact" # caption
  ) %>%
  kable_styling(latex_options = c("striped", "hold_position"), font_size = 10)%>%
  add_header_above(c(" " = 1, "Mean" = 3, "Difference" = 3, "p-value" = 3)) %>%
  cat(file="tab_uncontact_balance.tex")





# **Figure D1 Controlling for UN contact ----

# Control group
model<-lm(donation~fra+un+un_see+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg<-coeftest(model,vcov)
reg


new.data<-data.frame(group=c("C","UN","Fr"))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

# UN group
model<-lm(donation~untreat+un_see+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg<-coeftest(model,vcov)
reg

df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)


# **** France group 
model<-lm(donation~fratreat+un_see+factor(enumid)+factor(day), data=nomikos_jop_data)
summary(model)
vcov <- cluster.vcov(model, nomikos_jop_data$cluster)
reg<-coeftest(model,vcov)
reg

df<-model$df.residual
se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[3]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(reg[1,1]+est,0)

ggplot(aes(x=group, y=fit,ymin=ymin95, ymax=ymax95),data=new.data)+
  geom_linerange(size=1)+
  theme_bw()+
  geom_line(group=1,linetype=2)+
  ylab("Amount Sent to Tuareg Partner") + xlab("Treatment") +
  stat_summary(fun.data = mean_cl_normal, fun.args = list(conf.int = 0.95), geom = "errorbar",
               ymin=new.data$ymin95,ymax=new.data$ymax95, position = position_dodge(width = 0.90), width = 0.1) +
  geom_point(size=18,shape=21, fill=c("white", "lightblue", "grey"))+
  geom_text(x=c(1,2,3), y=new.data$fit, label=round(new.data$fit,2), size=5)+
  theme(text = element_text(size=14),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
        axis.ticks.x=element_blank(),
        axis.text.x = element_text(size=12))



dev.copy(pdf,"fig_app_rob_uncontact.pdf", height=6, width=8)
dev.off()



# **Figure D2: Women Only ----

model<-lm(donation~un+fra+factor(enumid)+factor(day), data=lab_female)
summary(model)
vcov <- cluster.vcov(model, lab_female$cluster)
reg<-coeftest(model,vcov)
reg

new.data<-data.frame(group=c("C","UN","Fr"))
new.data$group<-factor(c("Control", "UN","France"),levels=c("Control", "UN","France"))
df<-model$df.residual
se<-reg[1,2]
est<-reg[1,1]
new.data$ymin95<-est-qt(c(.975), df=df)*se
new.data$ymax95<-est+qt(c(.975), df=df)*se
new.data$fit<-round(est,0)

se<-reg[2,2]
est<-reg[2,1]
new.data$ymin95[2]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[2]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[2]<-round(reg[1,1]+est,0)

se<-reg[3,2]
est<-reg[3,1]
new.data$ymin95[3]<-reg[1,1]+est-qt(c(.975), df=df)*se
new.data$ymax95[3]<-reg[1,1]+est+qt(c(.975), df=df)*se
new.data$fit[3]<-round(reg[1,1]+est,0)

ggplot(aes(x=group, y=fit,ymin=ymin95, ymax=ymax95),data=new.data)+
  geom_linerange(size=1)+
  theme_bw()+
  geom_line(group=1,linetype=2)+
  ylab("Amount Sent to Tuareg Partner") + xlab("Treatment") +
  stat_summary(fun.data = mean_cl_normal, fun.args = list(conf.int = 0.95), geom = "errorbar",
               ymin=new.data$ymin95,ymax=new.data$ymax95, position = position_dodge(width = 0.90), width = 0.1) +
  geom_point(size=18,shape=21, fill=c("white", "lightblue", "grey"))+
  geom_text(x=c(1,2,3), y=new.data$fit, label=round(new.data$fit,2), size=5)+
  theme(text = element_text(size=14),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0)),
        axis.ticks.x=element_blank(),
        axis.text.x = element_text(size=12))



dev.copy(pdf,"fig_app_rob_women.pdf", height=6, width=8)
dev.off()
